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The multidimensional topography of the collective poten- 
tial energy function of a so-called strong glass former (silica) 
is analyzed by means of classical molecular dynamics calcula- 
tions. Features qualitatively similar to those of fragile glasses 
are recovered at high temperatures : in particular an intrin- 
sic characteristic temperature T c ~ 3500K is evidenced above 
which the system starts to investigate non-harmonic potential 
energy basins. It is shown that the anharmonicities are essen- 
tially characterized by a roughness appearing in the potential 
energy valleys explored by the system for temperatures above 
T e . 

PACS numbers: 61.43. Fs, 61.20.Ja, 64.70.Pf 

Even though the manufacturing of glasses by a quench 
from the high temperature liquid phase is a standard 
practice, the precise understanding of what happens at 
a molecular level in these materials, when lowering the 
temperature through the "glass transition" , remains an 
important theoretical challenge. Several approaches have 
emerged in the literature depending whether the struc- 
tural or the dynamical aspects are emphasized. In partic- 
ular the use of the concept of topological frustration [jlj||] 
or the use of the mode-coupling theory || proceeds from 
completely different points of view. A promising way to 
conciliate these approaches is the analysis of the shape of 
the potential energy function $(fi,f2, ■■■) in the multi- 
dimensional space of the coordinates of the interacting 
particles |||5|. It has been recently shown || that the 
topological analysis of the energy landscape of a fragile 
glass-forming liquid described by a two-body Lennard- 
Jones (LJ) potential allows to explain the distinct dy- 
namical regimes experimentally observed. 

In this letter, we present for the first time such an en- 
ergy landscape analysis in the case of a so-called strong 
glass former, namely vitreous silica, by means of clas- 
sical molecular-dynamics (MD) calculations. As in the 
case of fragile glasses, from a quantitative analysis of the 
inherent potential energy basins explored at various tem- 
peratures, one can define an intrinsic characteristic tem- 
perature, here equal to T c ~ 3500K, above which the 
system experiences anharmonicities. These findings are 
in agreement with a recent simulation study in which a 
change in the dynamical properties of a similar liquid sil- 
ica system has been observed around this temperature 
JtJ , suggesting that supercooled liquid silica behaves like 
a fragile glass former at temperatures above T c . Here we 



show that the anharmonicities exhibited by the energy 
basins explored at temperatures above T c are essentially 
due to a roughness in the shape of the basins, as previ- 
ously observed in simpler liquid systems Using an 
original procedure, we perform a quantitative analysis of 
this roughness and show that it is characterized by a typ- 
ical length (in the multidimensional configuration space) 
of about 0.2 A which is the signature of the annihilation 
of structural defects along the path down the potential 
energy valleys. 

Our silica system consists of 216 silicon and 432 oxygen 
atoms confined in a cubic box of edge length L — 21.48A 
which corresponds to a mass density very close to the ex- 
perimental value of vitreous silica. The constant energy, 
constant volume classical MD calculations have been per- 
formed using the sophisticated potential first introduced 
by van Beest et al. || and justified by ab initio calcula- 
tions. It has been recently shown to describe quite well 
structural |S,M, and vibrational jy], as well as relax- 
ational M and thermal jl2| properties of both, super- 
cooled viscous liquid and glassy silica. We have used the 
same input parameters and the same modus operandi as 
in our previous structural study After a full equili- 
bration of the liquid (about 28 ps) the system has been 
cooled down to zero temperature at a quench rate of 
2.3 x 10 14 Ks _1 which has been obtained by removing 
the corresponding amount of energy from the total en- 
ergy of the system at each iteration. At several temper- 
atures during the quenching process the configurations 
(positions and velocities) have been saved. With the use 
of these configurations as inputs of the MD calculations, 
the system has then been allowed to relax for a maximum 
duration of 84 ps after the quench (during this constant 
volume relaxation period the temperature of the system 
is only slightly increasing). The same procedure has been 
repeated for ten different initial liquid configurations and 
consequently all the physical quantities reported below 
always result from an average over these ten independent 
samples. 

For each collected temperature, immediately after the 
quench or after a given relaxation time (42 or 84 ps), 
the typical potential energy basin sampled by the sys- 
tem has been investigated using a procedure described 
by Delia Valle and Andersen Jp| which allows to perform 
a steepest descent from the initial multicomponent space 
configuration down to the closest underlying potential 
energy minimum, often called the "inherent structure" 
in the literature H . For that purpose a modified version 
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of the MD algorithm is adopted. At each MD step, the 
scalar product of the velocity by the force is calculated 
for each particle. If the product is positive the velocity 
of the particle is replaced by its projection on the force 
and otherwise the velocity is set to zero. During the de- 
scent, both the distance and the potential energy $ are 
calculated, allowing to determine precisely the shape of 
the energy basin. There are different ways to define the 
distance per particle in the multiparticle space. We have 
chosen to calculate either the direct distance x from the 
initial configuration, defined by |Q : 
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where the sum runs over all the particles located at Tt 
with a mass m il or the cumulated distance, X, calculated 
along the steepest descent path by : 
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where n labels the MD steps. The process is stopped 
when which is a monotonically decreasing function of 
x (resp. X), reaches a minimum $ m at x = x rn (resp. 
X = X m ). Practically we have chosen to stop the descent 
when - $("+!) < 1(T 6 eV. 

As a first result we present in Fig.l, the potential en- 
ergy minima <£>,„ obtained at different temperatures. The 
curve obtained immediately after the quench is similar 
to the one obtained previously with a different potential 
p3| . Of course, the system investigates energy basins of 
lower minima at lower temperatures. The slowing down 
of the decrease of $ m with decreasing temperature oc- 
curring between 4000K and 3000K is the signature of 
the glass transition as the system gets trapped after the 
quench in energy basins with almost the same minimum 
below a given temperature T g . This temperature is con- 
sistent with the estimate T g ~ 3500K obtained from a 
structural analysis done for the same system with the 
same quenching rate j^] but also with the extrapolated 
value obtained in a different study concerning the influ- 
ence of the quenching rate on the properties of the sys- 
tem [ [To| . Moreover, in agreement with the study done 
on the L J glass || , one observes in Fig. 1 that the inher- 
ent structure depends on the duration of the relaxation 
process after the quench. The relaxed curves exhibit a 
minimum around T g which is more and more marked as 
the aging time increases. This can be easily understood 
since for temperatures close to T g the system takes ad- 
vantage of the relaxation process in order to find lower 
energy minima, because it has enough kinetic energy to 
have a chance to cross the energy barriers between min- 
ima. Of course this chance becomes considerably smaller 
at lower temperatures in agreement with what is known 
from the thermal evolution of the relaxation time in the 



glassy phase of silica (note that in || the energy per par- 
ticle at T = OK depends on the cooling rate while we 
do not observe significant relaxation effects at T — OK). 
These results are consistent with a previous analysis of 
history effects in the same system |l5|] and also with the 
fact that T g should be smaller for a lower quenching rate 
]l0| . If we take into account both, the 15 ps necessary 
to reach 3500K from 7000K and the 84 ps of further re- 
laxation at 3500K, we obtain an effective quenching rate 
about six times smaller than the one we have used. Ac- 
cording to the dependence of T g with the quenching rate 
proposed by Vollmayr et al. |T0), this would correspond 
to a decrease of T g by more than 500K. 

In Fig. 2 we have plotted, as a function of temperature, 
x^ the square of the distance between a given initial con- 
figuration and the corresponding position of the inherent 
structure. This curve is remarkably similar to what has 
been previously obtained for Lennard- Jones glasses || 
and can be interpreted in the same manner. The squared 
distance x^ is linear with temperature up to a charac- 
teristic temperature T c ~ 3500K (here very close to T g ) 
above which anharmonicities appear. The same analysis 
as the one done in Q shows that T c corresponds to a 
change in the nature of the relaxation process, i.e. from 
diffusion above T c to hopping below T c . This is con- 
sistent with the recent work of Horbach and Kob on the 
same system Q who found a breakdown of the Arrhenius 
behavior of the transport coefficients above a tempera- 
ture close to 3300K and suggested that this characteris- 
tic temperature corresponds to the T c predicted by the 
mode coupling theory. A more striking result is that the 
curves in Fig. 2, in contrast with those in Fig.l, are not 
very sensitive to the duration of the relaxation process 
after the quench. This is also true for another charac- 
teristic quantity of the energy basins leading to the in- 
herent structures, namely A<!> = — $ m , the energy 
difference between the initial structure and the inherent 
structure, which has been plotted versus T in the inset 
of Fig. 2. It exhibits the same departure from the ex- 
pected harmonic linear regime (A$ = ^kgT represented 
by the dashed line in the inset) at T c and also the same 
remarkable stability against relaxation. This shows that 
as the aging time increases the system explores deeper 
and deeper basins (especially around T g ) as shown by 
the variation of $ m in Fig.l but the intrinsic characteris- 
tics of these basins, the mean width and the mean height 
(associated with x 2 ^ and <& < -°- ) — <f> m resp.), do not de- 
pend on their absolute position $ m in the energy scale. 
Therefore the characteristic temperature that we can de- 
fine here, T c ~ 3500K, is intrinsic and does not depend 
on the quenching rate. The fact that we obtain T g close 
to T c is simply due to our very large quenching rate. It is 
worth noticing that the same qualitative behavior than 
the one observed in Fig. 2 is also obtained by plotting X^ 
the square of the cumulated distance instead of x^- In 
that case the departure from a linear regime above T c is 
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even more pronounced as we find that the ratio X m /x m 
remains constant (of order f .2) for T < T c and increases 
markedly with temperature for T > T c . 

To investigate further the nature of the anharmonici- 
ties of the energy basins explored by the system for tem- 
peratures above T c , we propose a quantitative analysis 
of the curves F(X) = —d$/dX where X is the above 
defined cumulated distance (we have used — d<fr/dX in- 
stead of —d&/dx because it corresponds better to a local 
characteristic of the shape of the basins). Typical exam- 
ples of such curves are given in Fig. 3 for samples relaxed 
during 84 ps after the quench (as already shown in Fig. 2 
the aging time does not influence significantly the nu- 
merical results). It turns out that for T > T c , the $(X) 
curves, while always decreasing, exhibit step-like singu- 
larities evidenced by clearly visible peaks in the deriva- 
tive. Note that another manifestation of what can be 
called a roughness is the increase with temperature of the 
ratio X m /x m above T c , mentioned earlier (the same ob- 
servation has already been done for LJ glasses |J). To an- 
alyze quantitatively the roughness of the F(X) curves, we 
have first eliminated the overall mean general evolution 
by calculating the difference f{X) = F(X)- < F(X) > 
where < F(X) > is a local average of the data between 
X — SX and X + 5X (for convenience we have chosen 
SX = X rn /20 and limited the range of X values between 
4SX and X m — SX). The f(X) curves corresponding to 
the F(X) curves depicted in Fig. 3 are shown in the inset 
of the figure (they have been artificially shifted in the 
vertical direction for clarity). Subsequently the rough- 
ness of the curves f(X) has been analyzed by following 
a standard method [|l6| which consists in calculating the 
power spectrum S(k) defined as the Fourier transform of 
the autocorrelation product <?(£) = < f(X + £)f(X) > 
where the average is performed not only over the X val- 
ues but also over ten independent samples. The results of 
this analysis are reported in Fig. 4. In this figure we have 
reported as a function of temperature the mean intensity 
of the peaks measured by the standard deviation cry of 
the f(X) curves, which is equal to the square root of 5(0). 
Despite the error bars the curve exhibits a characteristic 
sigmoidal shape, indicating that the roughness only exists 
in the higher energy basins explored for T > T c and that 
the mean intensity of the peaks seems to saturate at very 
high temperatures. Furthermore, for T > T c , we observe 
that the power spectrum goes through a maximum and 
decays like k~ 1 after this maximum. This means that 
there exists a typical length (the inverse of the location 
of the maximum) X r characteristic of the mean distance 
between successive peaks while the curve between two 
peaks can practically be considered as "smooth". For 
T < T c there is no visible maximum in the power spec- 
trum which behaves roughly as k~ l over the whole k 
range. This is a further indication that the very weak 
roughness in the basins explored below T c has no signif- 
icance and that the low temperature basins can be con- 



sidered as smooth. The estimated values for X r above T c 
have been reported in Fig. 4. This typical length increases 
slightly with T from about 0.15 A at 4000K to about 
0.23 A at 7000K. It is interesting to relate the typical 
distance X r in the multidimensional configuration space 
to peculiar structural rearrangements occurring during 
the down-hill potential energy minimization process. On 
some specific high temperature samples we have com- 
pared the configurations between two successive peaks in 
F(X) and we have in each case observed the elimination 
of a single specific defect like a triconnected silicon atom 
or edge-sharing tetrahedra etc... In such rearrangements 
a "perturbed" cluster of about 30 to 50 connected atoms, 
is observed. Generally the largest displacement of about 
0.5 to 0.7 A is observed for an oxygen atom at (or very 
close to) the defect. Therefore the typical value of X r 
results from an average between the largest displacement 
near the defect and the "screening cloud" of displaced 
atoms connected to it. One can understand that X r be- 
comes insignificant below T c because the defects become 
rare in the low temperature basins as already shown in 

In conclusion we have numerically investigated the po- 
tential energy landscape of super cooled liquid silica de- 
scribed by the BKS potential using a steepest-descent 
molecular-dynamics scheme. We have shown that the in- 
herent structures sampled depend strongly on the effec- 
tive cooling-rate especially around J\ similarly to what 
was found in a Lennard- Jones glass [g| . Nevertheless at a 
given temperature the characteristics of the energy basins 
(mean height and mean width) seem to be insensitive to 
the history of the system. From a quantitative analysis 
of the potential energy valleys explored at various tem- 
peratures we have evidenced a characteristic temperature 
T c above which non- harmonic effects become dominant. 
This is consistent with an other recent study and in- 
dicates that strong and fragile glass formers are quite 
similar when studied near T c where a change in the na- 
ture of the relaxation process takes place in the liquid 
phase. We think that the distinct characteristics of a 
strong glass former appear mainly in the supercooled liq- 
uid and glassy phases below T c . Unfortunately the high 
value of T g that we obtain in our MD simulations does 
not allow us to study the temperature range between T g 
and T c . Furthermore using an original quantitative anal- 
ysis, we have shown that the anharmonic character of the 
higher energy valleys explored by the system above T c , 
is due to some roughness occurring in the shape of the 
potential energy basins. The existence of such roughness 
has already been invoked in the case of simpler systems 
HH but not quantitatively analyzed. In the case of sil- 
ica we have shown that this roughness is characterized by 
a typical length of about 0.2 A in the multi-dimensional 
configuration space, and can be associated with a sequen- 
tial elimination of defects when following the path leading 
down to the inherent structures. 
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Part of the numerical calculations were done at 
CNUSC (Centre National Universitaire Sud de Calcul), 
Montpellier. 
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FIG. 1. Plot of <J? m , the minimum value of the potential 
energy basin explored by the system, as a function of the 
temperature T down to which the system has been cooled 
from the liquid state. The three curves correspond to the 
different durations of the relaxation period after the quench: 
□ : no relaxation; o: 42ps; •: 84 ps. The data result from an 
average over ten independent samples. 
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FIG. 2. Plot of x 2 m , the mean square distance between the 
initial configuration and the corresponding minimum poten- 
tial location as a function of temperature. The three curves 
correspond to the different durations of the relaxation period 
after the quench: 

□ : no relaxation; o: 42ps; •: 84 ps. The data result from 
an average over ten independent samples. The straight line 
corresponds to a pure linear behavior with temperature. 
In the inset the variation of A$ = — <E> m , the energy 
difference between the initial configuration and the inher- 
ent structure, is plotted versus temperature with the same 
conventions (the dashed line represents the harmonic linear 
regime: A$ = |fe s T). 
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FIG. 3. 



Plot of the derivative F(X) = — ^ of some 
typical energy basins explored at different temperatures in 
the case of a sample relaxed 84 ps after the quench. From left 
to right the curves correspond to T = 1000, 3000, 5000 and 
7000 K, respectively. In inset are shown the corresponding 
curves f(X) obtained by making the difference between F(X) 
and a local average as explained in text. In the inset the 
curves have been arbitrarily shifted vertically for clarity. 
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FIG. 4. Plot of the mean-square intensity 07 (•) of the 
roughness in the space-derivative of the potential curve as a 
function of temperature after averaging over ten samples with 
a relaxation period of 84 ps after the quench. On the same 
figure the typical distance X r (□) between the peaks in f(X) 
is plotted as a function of temperature (right vertical axis). 
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